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Parampreet SinghB 

Department of Physics and Astronomy, Louisiana State University, Baton Rouge, Louisiana 70803, USA 

A brief review of various numerical techniques used in loop quantum cosmology and results 
is presented. These include the way extensive numerical simulations shed insights on the 
resolution of classical singularities, resulting in the key prediction of the bounce at the Planck 
scale in different models, and the numerical methods used to analyze the properties of the 
quantum difference operator and the von Neumann stability issues. Using the quantization 
of a massless scalar field in an isotropic spacetime as a template, an attempt is made to 
. highlight the complementarity of different methods to gain understanding of the new physics 

£N| ' emerging from the quantum theory. Open directions which need to be explored with more 

refined numerical methods are discussed. 

< 

^ ■ I. INTRODUCTION 

It is generally believed that a quantum theory of gravity will provide important insights on the 
resolution of singularities, the boundaries of spacetime in Einstein's theory of General Relativity 
(GR). Since we do not yet have a complete theory of quantum gravity, a pragmatic approach is to 
use known techniques of quantum gravity to understand quantization of spacetimes with reduced 
degrees of freedom, such as cosmological models. This provides a test bed not only for methods 
■ of full quantum gravity, but potentially opens new avenues to eventually link quantum gravity 

. with phenomenology and observations. An early attempt in this direction, based on Wheeler's 

ideas on geometrodynamics was Wheeler-DeWitt quantum cosmology. However, due to various 
1 mathematical and physical hurdles, in particular lack of guidance from a more complete theory, 

qq . it met with little success on the resolution of singularities. Though Wheeler-DeWitt quantization 

yielded consistent behavior in the infrared regime (i.e. at the classical scales) in agreement with GR, 
unfortunately its predictions agreed with GR even in the ultraviolet (large spacetime curvature) 
regime. As an example, a semi-classical state peaked in a macroscopic universe at late times in 
Wheeler-DeWitt theory, remains peaked on the classical trajectory throughout the evolution and 
^ I encounters the big bang singularity. 

In recent years, various shortcomings of Wheeler-DeWitt quantum cosmology have been over- 
come in the framework of loop quantum cosmology (LQC) [l|, which is a quantization of space- 
times with finite degrees of freedom based on loop quantum gravity (LQG), a candidate for non- 
perturbative canonical quantization of gravity. As in Wheeler-DeWitt theory, LQC is based on 
Dirac's method of constraint quantization. Instead of geometrodynamics variables in Wheeler- 
DeWitt theory, the classical gravitational phase space is labelled by the symmetry reduced versions 
of Ashtekar-Barbero connection A l a and its conjugate triad Ef, the phase space variables in LQG. 
Classical Hamiltonian constraint is expressed in terms of the elementary variables of quantization, 
the holonomies of connection and triads, and is quantized. An inner product is found using meth- 
ods of group averaging 0] which leads to a physical Hilbert space, and physics is extracted using 
a family of Dirac observables. This algorithm has been carefully carried out for various models 
in LQC. At the level of quantum Hamiltonian constraint, a striking difference between Wheeler- 
DeWitt quantum cosmology and LQC turns out to be the following. Unlike Wheeler-DeWitt 
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theory, where the quantum Hamiltonian constraint yields a differential equation on a continuum 
spacetime, the quantum evolution equation in LQC turns out to be a difference equation in geo- 
metrical variable(s). The origin of discrete quantum evolution operator lies in the discrete action 
of holonomy operators, which captures the non-local nature of field strength of the connection on 
the underlying quantum geometry. 

LQC started with the seminal works of Bojowald more than a decade ago, which indicated 
resolution of singularities, albeit at a kinematical level [1, Q]. These ideas were put on a firm 
mathematical footing by Ashtekar, Bojowald and Lewandowski in Ref. [jj. A loop quantization 
of cosmological spacetime, where the inner product, Dirac observables and the resulting physical 
implications using states in the physical Hilbert space were first made available, was performed for 
a spatially flat isotropic model with a massless scalar field [6H8J] . For this model, classical solutions 
are generically singular, and there is no possibility, as is sometimes true in models with potentials, to 
fine tune any parameter to avoid the singularity. Thus, it is an ideal arena to test whether methods 
of LQG can resolve the singularity. The massless scalar field also allows to overcome a conceptual 
difficulty to understand evolution in constraint quantization. In the absence of external time, scalar 
field serves the role of internal clock, and physics can be extracted using relational observables. 1 
A detailed analysis of resulting physics through sophisticated analytical and numerical methods, 
demonstrated that the big bang singularity is resolved and replaced by a quantum bounce when 
energy density (p) of the scalar field became equal to a critical value, p « 0.41/9pi anc k [8|. At small 
spacetime curvature (or large volume), the quantum Hamiltonian constraint, a difference operator 
with uniform discretization in volume, approximates the Wheeler-DeWitt quantum constraint, and 
hence agrees with GR. Robustness of results on genericity of bounce and maximum energy density 
were confirmed with an exactly soluble model in LQC [9[. Thus, quantum geometric effects in LQC, 
not only cure the singularity problem, but also lead to a physically consistent infrared behavior. In 
the early quantizations of LQC 0-0]) the latter was a serious problem, which was noted in different 
wor ks EBQ. As we will discuss later, this limitation is tied to the structure of quantum difference 
equation in the old quantization, which is uniformly discrete in eigenvalues of triads (proportional 
to area). It turns out to have difficulties with the von Neumann stability issues at large volumes 
and leads to a critical density which depends on the phase space variables in such a way that the 
bounce could occur even at very small spacetime curvature [7J]. 2 It turns out that for the isotropic 
models, the quantization which yields a difference equation with uniform discreteness in volume, 
which is also referred as the improved quantization, is the only one which results in the correct 
infrared behavior for various types of matter and a critical density which is a fundamental constant 
0- 

With the success of improved quantization of LQC [8], quantization of isotropic spacetimes with 
massless scalar field and positive spatial curvature 14, la], negative spatial curvature H, 17], and, 
positive Jis, LL9[] and negative cosmological constant [20] were performed. 3 Each of these models 



introduced non-trivial subtleties in the quantization procedure, and unlike the case of spatially flat 
model, they could not be solved exactly. Extensive numerical simulations in these models confirmed 
the robustness of results on bounce, 4 and revealed various new features of rich physics resulting 



1 In the absence of fields, this role can be played by geometrical degrees of freedom. An example is the case of 
vacuum Bianchi-I model, where predictions can be extracted by treating one of the scale factors as internal time. 

2 There are additional problems with such a critical density which will be discussed in Sec. V. 

3 In above works, numerical techniques to find physical states in the quantum theory were largely developed by 
Tomasz Pawlowski. 

4 In these works, occurrence of bounce was studied by computing expectation values of volume operator. Recently, 
bounce in massless scalar model has been established by computing quantum probabilities using consistent histories 
approach [21J. A similar analysis for Wheeler-DeWitt theory predicts probability of encountering a singularity as 
unity 
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from quantum geometry. Recently, a rigorous quantization of Bianchi-I model with a massless 
scalar field has been performed [ijj ]. which overcomes problems with an earlier quantization [24]. 
Quantization of Bianchi-II [2^] and Bianchi-IX models [2^] have also been proposed, and develop- 
ments on cosmological models, have also provided insights on loop quantization of Schwarzschild 
spacetimes 27H30I]. As in the case of isotropic models, to loop quantize these spacetimes rigor- 



ously and consistently requires overcoming several mathematical hurdles. From the perspective of 
numerical methods, one of the difficulties is that unlike the case of isotropic models, the quantum 
Hamiltonian constraint in anisotropic and black hole spacetimes turns out to be a partial differ- 
ence equation. With the non-uniform discretization in different variables, numerical analysis of 
loop quantum Hamiltonian constraints in these models is far more challenging than the isotropic 
models. 

In the above discussion we focused on the role of numerical methods following the analysis in 
Refs. 0-0]. However, use of numerical techniques to understa nd prop erties of quantum difference 
equation in old quantizations of LQC are dated earlier 31-37]. In these works, carried out 

before inner product and physical Hilbert space in LQC became available, the thrust of most of 
the analysis was to gain insights on viability of particular quantizations in LQC by using ideas 
of von Neumann stability to relate solutions of quantum difference equation at large volumes to 
the Wheeler- De Witt equation. These works also introduced notions of pre-classicality of solutions 
to identify the underlying conditions for a consistent infra-red behavior [L?l l32l ]. It is important 
to note some important differences in the usage of von Neumann stability analysis in LQC, and 
its conventional treatments in computational physics. Unlike the conventional treatments on von 
Neumann stability issues where the continuum partial differential equation (PDE) is fundamental 
and the goal is to understand whether a finite difference equation obtained by discretization of PDE 
provides a good approximation, in LQC it is the quantum difference equation which is fundamental 
and the goal of stability analysis is to verify whether in the large volume limit (or at small spacetime 
curvatures), solutions of quantum difference equation agree with those of the Wheeler-DeWitt 
equation. Note that the quantum difference equation obtained in LQC is not obtained by a 
discretization of the Wheeler-DeWitt theory, and the goal of stability analysis is not to compare its 
solutions with Wheeler-DeWitt equation at ultraviolet (or small volume) scales. Using these ideas, 
the application of von Neumann stability techniques has provided useful insights on the infrared 
behavior of quantizations of various models in LQC. Currently, techniques are being developed 
to overcome problem of non-uniform discretizations in partial quantum difference equations in 
anisotropic and black hole spacetimes in LQC 3JJ, |4fl|. 

Another avenue where numerical techniques in LQC have been extensively used is in the effec- 
tive dynamics approach. This is based on the observation that the underlying quantum dynamics 
in various models-^js approximated by an effective spacetime description derived from an effec- 
tive Hamiltonian 
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l42j|. In this approach, instead of using quantum difference equations, one 
works with the modified dynamical equations on a continuum spacetime encoding quantum gravity 
corrections. Thus, one can use standard numerical methods to solve differential equations. Never- 
theless, various studies in this approach require a tight control on errors which has been carefully 
achieved. Examples of such works include studies on dynamics in inflationary and Ekpyrotic po- 
tentials 43|, 44], computation of classical probabilities in inflation 45, 46], assisted and multi- field 
inflationary dynamics 47J , Gowdy models [48f] and various investigations in anisotropic spacetimes 
49-52]. Effective dynamics approach has provided important insights on various issues in LQC. 
These include effects on cosmological perturbations 53J], genericity of singularity resolution 5414561] . 
and constraints on possible discretizations of the quantum difference operator 13j, [57J, which shed 
light on issues noted in stability analysis of quantum difference equation. 

The goal of this manuscript is to provide a brief overview of various numerical methods used 
in LQC. These range from those yielding physical evolution in the quantum theory to those re- 
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lated to using von Neumann stability ideas on quantum difference equation in LQC. Due to space 
limitations, it is difficult to go in details of all the diverse models and their developments, and 
some omissions are unfortunately inevitable. These include detailed discussions on various results 
from effective dynamics where conventional numerical methods to solve differential equations suf- 
fice, factor ordering and symmetrization issues 5 , and in-depth discussions on pre-classicality issues 
which were developed before a rigorous formulation of LQC was available. Some of these results, 
such as issues of pre-classicality, are covered in detail in a recent review on numerical methods in 
LQC which provides a complementary view to various developments [59]. To make our discussion 
simpler to follow and to integrate results from different methods, we take the strategy of explaining 
different methods and results using a single model - quantization of spatially fiat isotropic model 
with a massless scalar, and provide summary of main results of other models. 

The plan of this manuscript is as follows. In Sec. II, we describe the basic setting of loop quan- 
tization using the example of a spatially flat Friedmann-Robertson- Walker model with a massless 
scalar field which was rigorously quantized in LQC in Refs. This model serves as an ex- 

cellent tool to illustrate various details of the quantization procedure and properties of quantum 
theory, and has served as a template for quantization of models with spatial curvature 1414171] , in 
presence of cosmological constant 18-2(3] and anisotropic models 23, 23, 2^, 60]. After describing 
the classical phase space in Sec. IIA, we briefly discuss the quantum theory in Sec. IIB and obtain 
quantum Hamiltonian constraint as a difference equation. We conclude this section with a brief 
discussion of Wheeler-DeWitt quantization and the classical limit of loop quantum evolution equa- 
tion. (Readers who are interested only in the numerical techniques can skip Sec. II and refer to 
eq.(2.18)). Sec. Ill is devoted to numerical methods used to obtain physical solutions of the quan- 
tum theory. We focus on the way numerical simulations are carried out in spatially flat isotropic 
model with a massless scalar field 0-0], and discuss two methods to obtain physical states: one by 
using a fast Fourier transform (FFT) and another using evolution in internal time in Sec. IIIA and 
Sec. IIIB respectively. We summarize the results of spatially flat isotropic model in Sec. IIIC, and 
the other models in Sec. HID. In Sec. IV, we discuss the von-Neumann stability analysis. After 
explaining various subtleties and differences with the conventional stability analysis, we discuss 
examples of massless scalar, positive cosmological constant in LQC, and difference equation in the 
old quantization in LQC [5j, |7J. The latter case serves to illustrate limitations of isotropic models 
which do not yield a quantum difference equation which is uniformly discrete in volume, such as 
the old quantization in LQC 0, 0] and lattice refined models which allow arbitrary discretizations 
[63 ] . Sec. V deals with a brief discussion of effective spacetime description of LQC, and the way 
effective dynamics sheds insights on issues in stability analysis and uniqueness of discretization in 
isotropic LQC. We conclude with a summary and outlook in Sec. VI. 



II. LOOP QUANTUM COSMOLOGY: SPATIALLY FLAT ISOTROPIC MODEL 

In this section we provide a brief overview of loop quantization of spatially flat Friedmann- 
Robertson- Walker (FRW) spacetime with a massless scalar field. Our discussion is based on the 
analysis in Refs. [6j-l9j . We start with a discussion of the classical phase space, relation of connection 
and triad variables with the metric variables, and the classical Hamiltonian constraint. We then 
discuss kinematical aspects of the quantum theory and demonstrate the way physical Hilbert space 
and Dirac observables are obtained. We obtain quantum difference equation which turns out to be 
uniformly discrete in volume, and discuss its classical limit, which agrees with the Wheeler-DeWitt 



5 Various results discussed in this manuscript are unaffected by these issues. For readers interested in a detailed 
discussion on these issues we refer to Refs. [Til. IH^]. 
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A. Classical framework 

The spacetime metric for the spatially flat {k = 0) isotropic homogeneous spacetime is given by 

ds 2 = -N 2 dt 2 + a 2 {t)d^ 2 (2.1) 

where N is the lapse function and a(t) denotes the scale factor of the universe. The spatial manifold 
for the k = model can be non-compact with M 3 topology, or compact with a 3-torus (T 3 ) topology. 
If the topology is non-compact, we must introduce a fiducial cell V as an infrared regulator to avoid 
divergences in spatial integrations in order to define a Hamiltonian framework. An obvious but 
important consistency requirement is that the choice of the cell V must not affect the physical 
predictions of quantities invariant under the change of the cell. 6 The fiducial volume of the cell V 
is given by V a = J v d 3 x y/q, where q denotes the determinant of the fiducial metric q ab on the spatial 
manifold defined by the co- moving coordinates x a . The physical metric on the spatial manifold q a b 
is related to the fiducial metric as q ao = a 2 q ao , and the physical volume of the cell V is given by 
V = a 3 K V3 . 

The gravitational phase space variables in LQG are the Ashtekar-Barbero SU(2) connection A a 
and its canonical conjugate, the densitized orthonormal triads Ef, which satisfy 

{4(a0,4(l/)} = 87rG 7 <W 3 (:r,y) . (2.2) 
Here 7 ~ 0.2375 is the Barbero-Immirzi parameter whose value is set by the black hole ther- 



modynamics in LQG [61(. The connection and the triad variables are related to the conven- 
tional metric variables in the following way. The triad Ef are related to the spatial metric as 
EfEj = qq ab , and the Ashtekar-Barbero connection is related to the extrinsic curvature K ao as 
A\ = r a + -fK l a = T i a + 1 e bi K ab , where r* denotes the spin connection which vanishes for the k = 
model, e ai is the undensitized triad, and K a is the extrinsic curvature 1-form. 

Given the symmetries of the isotropic and homogeneous FRW spacetime, the connection A\ 
and triad Ef can be expressed in terms of an isotropic connection and triad pair (c,p) fl: 

A\ = cV-Wul Ef = pV^yfqe? , (2.3) 

where ef denote fiducial triads and are the fiducial co-triads compatible with the fiducial metric 
qab. The canonical pair (c,p) satisfy the Poisson bracket relation {c,p} = 8irGj/3. Note that the 
triad can take positive or negative values depending on the relative orientation of the physical and 
fiducial triads. Since we will not consider fermions in our analysis, the choice of an orientation 
represents a gauge freedom which does not affect physical predictions. It will be later fixed by 
choosing symmetric states in the quantum theory. It is also useful to note the relation between 

(c,p) and the metric variables. The triad p is related to scale factor as \p\ = Vo^a 2 , and only on 

1/3 

the classical solutions of GR, c = 7V0 a/N, where a 'dot' denotes derivative with respect to time 
t. 



6 As will be discussed later, this requirement can be used as one of the criteria in selection of a unique discretization 
of quantum difference evolution operator in LQC [T3 |. 
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It turns out that quantum theory is considerable simpler if we express the gravitational phase 
space in terms of a canonical pair (b,u) related to (c,p) as, 



c Id 3 / 2 

where sgn(p) is +1 if physical and fiducial triads have same orientation, and is -1 if the orientation 
is opposite. The conjugate variables b and v satisfy, 

{ h M = \, (2-5) 

where h in the denominator is an artifact of the definition of v in (|2.4h . It should be noted that 
though v is a measure of physical volume V = \p\ 3 ^ 2 , it has dimensions of length. The variable b 
has dimensions of inverse length, and only in the classical theory is related to the Hubble rate H of 
the scale factor as b = jH = jd/a (for the lapse N = 1). The matter part of the phase space, for 
the massless scalar model under consideration, is labelled by the scalar field 4> and its conjugate 
momentum p^, which satisfy {4>,p<f>} = 1. The conjugate momentum is related to 4> as p^ = V Q a 3 cj) 
(for N = 1). The classical phase space is thus four dimensional identified by (6, v\ (j),p^). 

Before we proceed to the Hamiltonian constraint, it is useful to note the way these variables 
transform under freedoms related to the choice of spatial coordinates and fiducial volume V Q which 
brings forward the advantage of using b instead of c in the framework. The first freedom allows 
a rescaling of coordinates x — > x' = ax without affecting the metric. Under this freedom, various 
variables considered above transform as 

a -> a~ l a, V D ->■ a 3 V Q , V V, (c,p) -> (c,p), (b,u) ->■ (b,u) and -> (4>,p^) ■ (2.6) 

Thus, under the rescaling of coordinates, both c as well as b are invariant. However, under the 
freedom of the choice of fiducial cell, V — > /3 3 V, one obtains 

a ^ a, V -> pVo, V -> /3 3 V, (c,p) -> (/3c,/3 2 p), (6,^) -> (6,/3 3 ^) and (<P, P4> ) -> (0,/3%) . 

(2-7) 

It is to be noted that under the rescaling of fiducial cell, b is invariant but c is not. Further, it is 
straightforward to check that if we consider a phase space variable P m = c\p\ m , then it is invariant 
under rescaling of fiducial cell only when m = —1/2, i.e. when P m = b. The variable P m is used 
in the lattice refinement scheme with —1/2 < m < [62] 7 , and this observation will be useful for 
later discussion (in Sec. V) to gain insights on the consistency of improved quantization [8j and 
inconsistencies of other quantizations in isotropic LQC. 

Before we write the classical Hamiltonian constraint in terms of gravitational phase space vari- 
ables, we note that the scalar field satisfies □</> = 0, and is used as internal time in LQC. It is 
then natural to choose the lapse N = a 3 such that the time variable r satisfies Dr = 0. Unless 
noted specifically, in the following, we will work with this choice of lapse function. Using this lapse 
function, the classical Hamiltonian constraint for the isotropic homogeneous model with a massless 
scalar field becomes 8 



The usage of the term 'lattice' here should not be confused with the one in lattice loop quantum cosmology in Ref. 
[63( | . where it is used in the context of a lattice of homogeneous universes to investigate inhomogeneities. 
The presence of h in the following classical constraint is an artifact of the definition of v in eq. (|2.4p . 
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Using the classical Hamilton's equations, it is straightforward to obtain the dynamical equation 
from the above Hamiltonian constraint, 

= ±Vl27rGV (2.9) 
whose integration yields the classical trajectories: 

' In — + (f> c (2.10) 



where u c and 4> c are constants of integration. In the classical theory, we thus obtain expanding 
and contracting branches. The expanding branch encounters big bang singularity at v = in the 
backward evolution at a certain value of internal time 4>. Similarly, the contracting branch ends in 
big crunch singularity in the future evolution. 



B. Quantum framework 

In LQG, elementary variables are given by the holonomies of connection A l a and fluxes of triad 
Ef . In LQC, to define holonomies we consider straight edges fie^. The holonomy of a connection 
c turns out to be, 

„<<" = co.(f)l-* sin (2.11) 

where I is a unit 2x2 matrix and Oi denote Pauli matrices. Similarly, fluxes are computed across 
the face of the cell V. Due to homogeneity, the flux integral turns out to be proportional to triads 
0]. Thus, the elementary variables for quantization are the elements of holonomies Nu := e l ^ C//2 
and triads p. Since fi is arbitrary, Nu generate an algebra of almost periodic functions of c. The 
gravitational part of the kinematical Hilbert space %kin = ^fjjT ® ^Sn** * s obtained by finding a 
representation of holonomy flux algebra using Gel'fand-Naimark-Segal construction, which results 
in %|[q V as a space of square integrable functions on the Bohr compactification of the real line: 
L 2 (MBohr> d/iBohr) M- Unlike the gravitational sector, the matter part of kinematical Hilbert space 
is obtained by quantizing matter with methods of Fock quantization. In contrast to the Wheeler- 
DeWitt quantization, where the gravitational part of the kinematical Hilbert space is L 2 (M, dc), 
normalizable states in 'H^ 1 ^ / can be expressed in terms of a countable sum of N^, with (Np\N') = 
6„„r . It should also be noted that unlike Wheeler-DeWitt theory, there is no operator corresponding 
to c (or its curvature) in LQC. Instead, its information is contained in the corresponding holonomy 
operator. The action of the operators on states ^(/u) in the triad representation is translational: 
N^fy(fj,) = ^(^u + £), where £ is a constant. In contrast, the action of c in Wheeler-DeWitt 
theory is differential in the triad representation (as in the standard Schrodinger mechanics). Thus, 
even at the kinematical level, there are significant differences between LQC and Wheeler-DeWitt 
quantizations. The most important difference being that the continuum differential geometry 
common to GR and Wheeler-DeWitt theory is replaced by a discrete quantum geometry in LQC. 

The next step in the quantization procedure is to express the classical Hamiltonian constraint 
in terms of the elementary variables for quantization. In terms of triads and holonomies, the 
gravitational part is given by 

C grav = 7 - 2 V- 1 / 3 e* jk el e) \p\ 2 Fj . (2.12) 
Here F a ^ is field strength of the connection, expressed in terms of holonomies over a square pla- 

1/3 

quette with length ^lV 

Fj = -2 lim Tr ~ 1 T A SAM , (2.13) 

ab ArD-»o V ArD / a b K 1 
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where h Ur] = hf ] hf \hf )-*. Before one promotes C grav to a quantum operator, it is 
important to note that the limit □« — > does not exist in the quantum theory. The non-existence 
of this limit is tied to the underlying quantum discreteness in LQG, and is consistent with the fact 
that in LQG there exist no operators corresponding to connection or its curvature. The underlying 
quantum discreteness allows the loop to be shrunk only to a minimum area, given by the lowest 
non-zero eigenvalues of the area operator in LQG, Alp\ where A = 4\/3irj. This constrains the 
parameter p to j23j 

A/ 2 

f = =f • (2-14) 
\P\ 

The dependence of p, on triad p makes the action of resulting field strength operator, 

Fjnm = ^v - 2/ *m{^^) m . (2-i5) 

less straightforward than the simple translation action of elements exp(zcrc/2) on the eigenstates 
of the triad. It turns out that the action simplifies in the volume representation (y). In the 
eigenbasis of v operator, the action of elementary variables V and its conjugate exp(iA6), where A 
is a parameter with dimensions of length, is given by 

V\u) = 2-!rylp l \i/\\v), exp^Xb)\v) = \v - 2A) . (2.16) 

The parameter A is related to the area gap as A 2 = A/|j. It thus directly captures the discreteness 
of the underlying quantum geometry. 

The resulting action of the gravitational part of the quantum Hamiltonian constraint yields a 
difference equation in uniform steps of volume. For the case of a massless scalar field, the action 
of total Hamiltonian constraint operator Ch = C gTSCV + 167rGC ma tt is given by, 

dpM) = frGu^v^VM) =: -8*(i/,0 (2.17) 

where O is a positive definite, internal time independent difference operator, with the following 
action: 

ev(v, <j>) := v ((v + 2A)^(^ + 4A) - 2z^(^, 6) + (u - 2X)^(u - AX)) . (2.18) 

4A Z 

The physical states ^f(i/,<p) are chosen to be symmetric under the change in orientation of triads, 
i.e. ty{y, (j)) = ty(—v,<f>). In the absence of fermions, as is the case in this analysis, this symmetry 
requirement eliminates a large gauge freedom induced by the parity operator II with the action 
n^I'(^) = u), and associated with the choice of the triad orientation. Note that the form of 
the quantum constraint is very similar to as in the Klein-Gordon theory, where 4> plays the role 
of time and acts like a spatial Laplacian operator, and as in the Klein-Gordon theory, physical 
states can be decomposed into positive and negative frequency subspaces. 

In order to extract physics, we need to find an inner product. It can be obtained by using group 
averaging procedure [2|, or demanding that the action of Dirac observables, V\$, the volume at 
internal time (p, and be self-adjoint. The physical inner product turns out to be 

(*i|* 2 ) = ^lfo&OM -1 *^^) • ( 2 - 19 ) 
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The action of Dirac observables on states ^f(u,(f)) is given by 

V\^{v,4>) = 27r 7 Z| 1 e i ^-^)| I /|*(^^ ), and p^(u,<f>) = -ihd^ 4>) = Vetf(v, <f>) . 

(2.20) 

The action of the volume observable V\^ can be understood as considering the state at <fi = <f) , 
multiplying by volume V, and evolving the state to 4>. Dirac observables preserve the positive 
and negative frequency subspaces, and thus it suffices to consider positive frequency states as the 
physical states, which satisfy —ihd^iy^cj)) = v9$(i',^). 

This provides us the physical Hilbert space 7i p h ys which is the space of normalized positive 
frequency states ty(v, 4>) which are symmetric under the change in orientation of triads. Since O 
is a difference operator, physical states have support on lattice v = dbe + 4nA, where e S [0,4), 
and the Hilbert space is decomposed into subspaces % e which are preserved by the evolution. 
To extract physics we can restrict ourselves to any of the sectors e and the results turn out to 
be independent of this choice. From the point of view of the big bang singularity, the most 
interesting sector is e = 0. (Numerical results on singularity resolution discussed in next section, 
correspond to this sector). With the availability of a self-adjoint quantum evolution operator, 
Dirac observables, and the inner product, we have a rigorous framework to extract physics from 
this quantum theory. 

The exercise carried out above can be repeated in a straightforward way for the Wheeler- 
DeWitt theory 0-0] • One can find the kinematical Hilbert space, action of quantum Hamiltonian 
constraint and Dirac observables V\$ and p^, and the inner product. A crucial difference is that 
unlike in LQC, the underlying geometry in Wheeler-DeWitt theory is not discrete and the quantum 
constraint (analogous to (|2.17|) , turns out to be @, 0] 

djply, <f>) = -9 <j>) := 12^G v d v v d v V(u, <f>) . (2.21) 

Thus, unlike LQC, where the evolution equation is a difference equation with uniform steps in 
volume, in Wheeler-DeWitt theory the evolution equation turns out to be a differential equation. 

We now show that the LQC quantum difference equation leads to the Wheeler-DeWitt equation 
in the large volume limit. For smooth wavefunctions, at v 3> lp\, we can expand ^{y ± 2A, (f>) in 
eq. flZHED as, 

± 2A, 0) = 0) ± 4A <f>) + \ (AX) 2 d?V(v, <j>) ±\ (4A) 3 %V(u, cf>) + ... (2.22) 

2 b 

Substituting in eq. (|2.17p . we obtain, 

3\ <P) ~ 12vrG v {Br&iy, <P) + ud^(u, 0)) + O (X 3 d 3 u V(v, 0)) . (2.23) 

If the wavefunction is slowly varying, in the approximation that the terms of the order A 3 ^^ can 
be neglected, we obtain eq. (l2.2ip . Thus, in the limit of large volume (small spacetime curvature), 
the quantum Hamiltonian constraint in LQC which is a quantum difference equation, yields the 
Wheeler-DeWitt equation. 

III. NUMERICAL TECHNIQUES FOR OBTAINING PHYSICAL STATES 

In the previous section, we showed the way loop quantization of a spatially flat isotropic model 
results in a quantum Hamiltonian constraint which is a difference equation with uniform discrete- 
ness in volume. Here we discuss numerical methods to determine physical states satisfying the 



9 



quantum Hamiltonian constraint [7j, l8[ . We first describe a method based on calculation of sym- 
metric eigenfunctions and computing a fast Fourier transform (FFT) to determine physical states. 
This is followed by description of another method, involving specification of an initial state peaked 
on a classical trajectory at late times in a macroscopic universe and evolved in internal time <f> us- 
ing quantum evolution equation (|2.17p . We summarize results from these methods, and for other 
models in Sec. IIIC and Sec. HID respectively. 

A. Physical states by FFT method 

This method is based on obtaining physical states by noting that any symmetric state in the 

(s) 

physical Hilbert space can be expressed in terms of symmetric eigenfunctions el of the operator: 

= dk^(k)e[ s) , (3.1) 

J — oo 

where ty(k) is a suitable function to define state profile, such as a Gaussian function. By numeri- 
cally finding symmetric eigenfunctions, one can perform the FFT in the above equation and obtain 
the physical state. Since this method relies on the properties of eigenfunctions, we discuss some of 
the features of the latter below. Symmetric eigenfunctions are obtained by considering symmetric 
combinations of eigenfunctions of the operator: 0e w (z^) = ui (y), where uj 2 = UirGk 2 . The 
eigenfunctions e u {y) have a 2-fold degeneracy on the Hilbert space %f5 av (or on Hf^T), a sub- 
space of physical Hilbert space with states having support on the lattice C\ e \ (or £_i e i). The way 
symmetric eigenfunctions are constructed depends on the choice of the sector (e). For e = 0or2, 
lattices C\ e \ and C\ £ \ coincide with each other under the transformation: v — > —u, and symmetric 
eigenfunctions can be constructed from either of the lattices. Whereas, for £^0,2, these lattices 
are distinct, and hence symmetric eigenfunctions are constructed using eigenfunctions e w {v) on 
positive as well negative lattices. In this case, eigenfunctions have a 4-fold degeneracy and sym- 
metric eigenfunctions are constructed by a linear combination of e-uju on positive and negative 
lattices. 

To obtain these symmetric eigenfunctions, we first compute eigenfunctions e u {v) numerically 
for a given e. Since at large volumes (V 3> lp{), 0—^0, eigenfunctions e w {v) are expected to 
approach the Wheeler-DeWitt eigenfunctions (^(v). This indeed turns out to be true. In Fig. 1, 
we show the behavior of an eigenfunction in LQC and compare it with that of the Wheeler-DeWitt 
theory. Here we note that each eigenfunction §^{v) of the Wheeler-DeWitt equation (|2.2ip can be 
expressed in a linear combination of basis functions e\ k \(is) as, 

&» = Ae\ k] (v) + Be_ ]k] {v) , (3.2) 

where A and B are constant coefficients, and at large volumes e u {u) ~ §^iy) + 0(\ 2 /v 2 ). It is 
interesting to note that the eigenfunctions e^{y) decay exponentially at small volumes. The volume 
at which exponential decay occurs is numerically found to be proportional to the value of uj (or 
the choice of p^). Recently, the exponential damping of eigenfunctions in the quantum regime has 
been studied analytically, which provides additional insights on the behavior of eigenfunctions {64j. 
It turns out that the volume at which exponential decay occurs corresponds to the volume at the 
bounce of the universe in LQC. 

After obtaining symmetric eigenfunctions, we need to choose a function ^(k) in the integral 
(|3.ip . We are interested in those states which at late times (or large volumes) correspond to a 
universe with a small spacetime curvature described by GR. We consider a sharply peaked state in 
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FIG. 1: An eigenfunction of quantum Hamiltonian constraint in LQC is compared with that of the Wheeler- 
DeWitt theory (solid line), for uj = 50. Here v corresponds to v/\ in Planck units. Eigenfunctions approach 
each other at large volumes. At small volumes, eigenfunction in LQC decays exponentially. 

volume as well as its conjugate b at late times on a classical trajectory. Due to relation between b, 
P0 and V implied by the classical constraint (|2.8p . it suffices to choose a state which is peaked in 
volume and p^. To be concrete let us choose such a state at a large value of p^ in the expanding 
branch: \p^\ 3> yGft. If k* (which is negative in the expanding branch) corresponds to the value 

of p^, via P0 = —Vl^nGhk, at which the state is peaked then for a small dispersion a we choose 

$(jfe) = e -{k~k'fl2a* e -iu<f>* _ (3.3) 

The volume at which state is peaked at time (ft = (ft Q is determined by the value of (ft* on using the 
equation for classical trajectory (|2.10p . To compute ^(^, (ft), one finds symmetric eigenfunctions in a 
range, say k* — Wa, k* + IOcj, by discretizing the interval in a large number of points (approximately 
2 12 ). The physical state ^{v, (ft) is computed by a FFT, and, the norm ||\P|| 2 and expectation values 
of Dirac observables, V\<j> and p<f> can be evaluated. We discuss the resulting physics in Sec. IIIC. 



B. Physical states by evolution in internal time 

This method deals with the specification of an initial state peaked at a large values of p§ 
and volume at (ft = cp with suitable boundary conditions, and evolve it in the internal time (ft 
using the quantum difference equation (|2.17p . 9 In this approach, a numerical computation of 
symmetric eigenfunctions is not required. Rather one has to numerically solve a large number of 
coupled equations in (ft, in a finite domain of numerical integration. To specify boundary conditions, 
one considers the quantum constraint equation in the form id^(u,(ft) = yQ^S(u,4>), which at 



For a detailed discussion of analytical issues involved in this construction, we refer the reader to Refs. M, |65| . 
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FIG. 2: The plot shows a section of the evolution of the amplitude of a physical state \E'(^, 4>) obtained using 
internal time evolution described in Sec. IIIB, and the resulting expectation values of the volume observable. 
The initial state is a solution of the Wheeler-DeWitt equation with a Gaussian profile peaked at v* — 22400 
and p*^ = 1500 (in Planck units), and is evolved using difference equation in LQC. Unlike the result of the 
state evolved in Wheeler-DeWitt theory (shown in Fig. 3), the state in LQC does not encounter singularity, 
but bounces when energy density becomes equal to p max ~ 0.41/jpianck. The bounce in this simulation occurs 
at vj, « 488, where Vf, denotes value of v at which bounce occurs. At small spacetime curvature, expectation 
values of volume observable are peaked on disjoint classical trajectories shown by solid lines. 



large volumes is approximated by d^(y,<f)) = Vl2TrGud u ^/(u, <f>). Using a discretized version of 
this equation, boundary conditions are specified permitting only outgoing solutions. In all the 
numerical simulations carried out in Refs. 0,0], the boundary was chosen such that the value of 
the wavefunction at the boundary was negligible compared to its peak value. 

In numerical simulations of spatially flat isotropic model, three methods have been used to 
specify initial data ^|^ o and d^\ ( f lo at (p = <j) . These are 0,11]: 

1. A Gaussian state normalized with respect to inner product in the Wheeler-DeWitt theory, 
which is peaked in p^ and volume for a macroscopic universe. The initial state is a minimum 
uncertainty state in V and b. 

2. A solution of Wheeler-DeWitt equation is chosen as an initial state with Gaussian profile 
peaked at a classical trajectory at small spacetime curvature. It is chosen to minimize 
uncertainty in 4>,p^. 

3. A similar method as the previous one with a variation taking into account the fact that 
eigenfunctions of in large volume limit satisfy eq. ()3.2p . To obtain this behavior, initial 
state is constructed by a phase rotation of eigenfunctions of the Wheeler-DeWitt theory. 
The phase required for a suitable rotation is determined numerically. 

With the above initial data, numerical simulations were performed for a large range of 
ranging from 100-10000 (in Planck units), for varying dispersions (~ 1% - 10%) in Refs. 0,Sj. 
More recently, these simulations have been carried out for a much smaller values of p^ (~ 5 - 20) , 
and higher dispersions (~ 15% - 30%) [6^]. In the numerical simulations performed in Refs. 0,0], 
numerical errors in internal time discretization were controlled by refining 'time' step such that 

ll^i - #2 1| < ||*||2(A0)e (3.4) 
where e is a very small number, and, and \&2 correspond to wavefunctions computed with 



two discretization of 4>: A(p and A(p/2. In simulations carried out in Ref. 66f] , consistency of the 
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FIG. 3: The evolution of initial state chosen in Fig. 2 is shown in Wheeler-DeWitt theory for comparison. 
The state peaked on the classical expanding trajectory evolves to the big bang singularity and does not 
bounce. The expectation values of volume observable are also shown for two different simulations, one 
corresponding to expanding branch and the other to the contracting branch. The states in expanding and 
contracting branches in Wheeler-DeWitt theory are peaked on classical trajectories (shown by solid lines), 
which are singular and disjoint. 



discrete time step has been found by taking into account Courant-Friedrichs-Lewy (CFL) condition 
for given properties of the initial state. In various simulations in different works, a typical evo- 
lution between <f> = to (j> = —2 is typically performed by discretizing time interval in 10 5 -10 6 steps. 

We conclude the discussion of Sec. IIIA and Sec. IIIB with a brief comment on the two methods 
to obtain physical states. Though both methods can be equivalently used to obtain physical states, 
they may have some advantages and disadvantages depending on the model considered. As we 
see above, unlike the method to find physical states by FFT discussed in Sec. IIIA, one needs to 
carefully address the issues of domain of dependence, convergence and stability in obtaining physical 
states by the method evolution. Thus, in comparison to the FFT method where physical states are 
obtained once one has obtained eigenfunctions, this method may appear more involved. However, 
a notable advantage of evolution method lies in the fact that it is computationally economical since 
it has less memory requirements in comparison to the FFT method. A naive implementation of the 
FFT method for spacetimes with more degrees of freedom, such as Bianchi-I model, can raise the 
computational cost significantly. On the other hand, the latter method, at least for some simple 
models, is comparatively straightforward to be parallelized. 



C. Key features of results and their robustness 



We summarize results from the numerical methods described above to obtain physical states 
in spatially flat isotropic model with a massless scalar. We also discuss, in brief, results from the 
spatially curved models with massless scalar [3, [H, 17], with positive and negative cosmological 



constants 



19 . l20l ] and from preliminary works in the presence of anisotropics. 



• The main result obtained using numerical simulations based on FFT and state evolution 
methods, is that physical states in LQC do not encounter classical big bang singularity. 
Rather, expectation values of the volume observable, {V\^} reveal that universe bounces from 
a forward evolution in contracting branch to an expanding branch (or vice-versa) when the 
energy density of the scalar field becomes p 0.41ppi anc k. The bounce occurs independent 
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of the choice of p$ or the sector e. This behavior can be seen in Fig. 2 where amplitude of 
the wavefunction in LQC is plotted versus v = v/\ and 4> and expectation values of volume 
observable are compared with classical trajectories. The initial state is peaked at v* = 22400 
(where v = v/X) with p*^ = 1500 (in Planck units), and evolved using method to consider 
initial state as a solution of Wheeler-DeWitt equation in Sec. IIIB. The lattice in this 
simulation is chosen to be e = 0. To compare with the Wheeler-DeWitt theory, we show the 
evolution of the same initial state specified in the expanding branch using Wheeler-DeWitt 
equation (|2.2ip in Fig. 3. Unlike the state in LQC which bounces, one finds that the state in 
Wheeler-DeWitt theory is peaked on the classical trajectory at all times and encounters the 
big bang singularity at v = in the past evolution. Note that unlike, Wheeler-DeWitt theory 
where expectation values agree with classical solution at all times, in LQC expectation values 
show a non-singular bounce which joins two disjoint trajectories of classical GR. 

• Classical GR turns out to be a good approximation to loop quantum evolution till the energy 
density reaches a percent of the Planck value. Significant deviations from classical trajectory 
occur when energy density becomes larger than this value. In regime where energy density 
is smaller than above value, a physical state in LQC is peaked on a classical trajectory in 
the expanding branch after the bounce, and on a classical contracting branch in pre-bounce 
regime, corresponding to the same value of p^. Quantum gravitational effects bridge two 
disjoint classical solutions. Thus, LQC yields GR in the infrared limit and new Planck scale 
physics devoid of singularity in the ultraviolet regime. This behavior can be seen from plot 
of expectation values in Fig. 2, which shows that at large volumes, there is an excellent 
agreement between LQC and classical GR, where as at small volumes there are significant 
departures. 

• Physical states remain sharply peaked through out the evolution, as is also clear from Fig. 2. 
Depending on the method used and properties of initial state, there is only a little variation in 
behavior of relative fluctuations across the bounce. Subsequent to the analysis in Refs. @,@], 
strong constraints on change in relative fluctuation of volume and pj, observables across the 
bounce have been obtained analytically for states which are similar as well as more general 



than those considered in numerical simulations 67H69I] . These analytical results confirm the 
numerical observation that if one considers an initial state peaked in a classical macroscopic 
universe at late times, it evolves via LQC quantum constraint to a state across the bounce 
with similar features at early times. 

• Norm of the state is preserved and expectation value of p^ observable remains constant 
throughout the bounce. 

• The expectation value of volume observable at the bounce coincides with the value at which 
eigenfunctions become exponentially damped 0, 0]. This result has also been confirmed 



analytically 64] . 



• Extensive numerical simulations confirm that underlying quantum dynamics in LQC can 
be approximated by an effective spacetime dynamics even in the Planck regime for states 
which peak in a macroscopic universe at late times. This behavior is shown in Fig. 4 and is 
discussed in detail in Sec. V. 

Remark 1: Spatially flat isotropic model with a massless scalar field can be solved exactly in the b 
representation 0] . Results obtained analytically confirm with those from the numerical simulations 
in the v representation performed earlier. In particular, using exactly soluble model one finds that 
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the bounce occurs for generic states in the physical Hilbert space, not necessarily semi-classical 
ones, in the sense that the expectation values (V\^,) are bounded below. The energy density 
of the scalar field has a finite upper bound given by p m ax = 

3/(8ttG 7 2 A 2 ) « 0.41/Opianck, which 

agrees with the value of critical energy density at which bounce occurs in all numerical simulations. 

Remark 2: In the old quantization of LQC 0, [2^] , the quantum difference equation turns out to 
be uniformly discrete in the eigenvalues of triads. This occurs if while relating area of the loop 
over which holonomies are considered to the minimum eigenvalue of area operator, p, is considered 
as a constant (// ). The resulting action of the elements of the holonomies exp(i^ c) is by uniform 
translations on triad eigenvectors, a feature shared by the action of the quantum Hamiltonian 
constraint in the old quantization. A rigorous quantization of this model was performed in Ref. 
@,0|. As in the improved quantization, there is a bounce, but with an important difference. Unlike 
the bounce discussed above, in the old quantization of LQC, bounce does not occur at a universal 
value of energy density. For the massless scalar model, the energy density at the bounce turns out 
to be inversely proportional to p^. Also, recall that p^ is not invariant under rescaling of fiducial 
cell (eq. (|2.7p ). 10 These features lead to a serious infrared problem in this quantization. The same 
problem plagues loop quantizations of isotropic models which do not yield quantum difference 
equation with uniform discreteness in volume (such as lattice refinement scheme [62]). We discuss 
these issues further in Sec. V. 

D. Other models 

In this part, we briefly summarize various numerical results obtained in spatially curved 
model, spatially flat model with positive and negative cosmological constant, and in presence of 
inflationary potential, and Bianchi-I model. 

Spatially curved models: Following the quantization procedure outlined in Sec. II and numerical 
techniques discussed above, a rigorous anal ysis of spatially closed (k = 1) isotropic model with a 



massless scalar field was performed in Refs. 14l.ll5l|. Unlike the k = model, a classical k = 1 filled 
with massless scalar inevitably collapses in future evolution of the expanding branch and encounters 
a big crunch singularity in a finite time. Thus, there is a physical singularity both in the past and 
the future evolution. LQC resolves both of these singularities for arbitrary choices of parameters. 
A technical difference in comparison to the k = model is that the corresponding quantum 
difference operator, ®k=i (which is uniformly discrete in volume) has discrete eigenvalues, due to 
which numerical procedure is more involved. Results from the numerical simulations with states 
peaked at classical trajectories at late times and evolved with the quantum difference equation 
show bounce at p « 0.41ppi anc k before big bang and big crunch are reached. In this model, loop 
quantum dynamics exhibits a cyclic evolution. Notably, states remain sharply peaked throughout 
the evolution for many cycles [[14]. As in the case of spatially flat model, an excellent agreement of 
loop quantum dynamics and effective dynamics obtained from an effective Hamiltonian was found 



1|. 



Loop quantization of k = 1 model overcame various limitations of an earlier quantization of the 
same model in LQC 0], which led to a quantum difference equation with uniform discreteness in 
triad eigenvalues. These limitations were first noted in the work of Green and Unruh [hJ, who 
found that eigenfunctions of the quantum constraint operator do not exhibit consistent infrared 



10 As discussed in Sec. V, the infrared problem occurs in this quantization irrespective of problem of rescaling of the 
fiducial cell. 
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properties. In particular, they do not decay at large volumes where recollapse of the classical 
k = 1 universe is expected. These problems were resolved in the quantization proposed in Ref. 



l4 . 1 151 ] . We will later discuss the way this issue is related to von Neumann stability of the quantum 



difference equation in the old and improved quantizations in LQC. 

The strategy for loop quantization of k = and k = 1 models runs into some technical 
difficulties for k = — 1 model [if], 03] • The existing quantization of this model considers holonomies 
of extrinsic curvature rather than the connection, and the resulting quantum difference operator 
is not essentially self-adjoint. Nevertheless, physical states have been studied by computing the 
eigenfunctions and taking a FFT (l6| . These states reveal existence of bounce when energy density 
becomes equal to p max ~ 0.41/9pi anc k. The quantization also shares various nice features of k = 
and k = 1 models, including states remaining sharply peaked at all values of internal time, and 
consistent ultraviolet and infrared behavior. Effective dynamics also provides a consistent picture 
of the evolution of physical states in LQC. 

Cosmological constant: Spatially flat model with a massless scalar field has been loop quantized 



with a positive 19|] and a negative cosmological constant 20]. As in the spatially flat case with a 
massless scalar, the difference equation turns out to be uniformly discrete in volume. For the case 
of negative cosmological constant, classical evolution predicts a recollapse of the universe at large 
volume, similar to the case of k = 1 model, and apart from past big bang, there is a big crunch 
singularity in the future evolution. The quantum difference operator turns out be self-adjoint, 
with discrete eigenvalues and its spectrum is numerically found to be non-degenerate. Numerical 
simulations reveal existence of bounce when total energy density of scalar field and cosmological 
constant become equal to p « 0.41/9pi anc k, and the big bang and big crunch singularities of the 
classical theory are avoided. Numerical studies with a large range of cosmological constant confirm 
that the recollapse of the universe occurs at the classically predicted volume. Further, states are 
found to be sharply peaked throughout the evolution, and effective dynamics is found to be in 
excellent agreement with quantum evolution [2o| . 

A rigorous quantization of spatially flat model with massless scalar and positive cosmological 
constant reveals various interesting features of the quantum difference operator [3, . In the 
classical theory, for any choice of positive cosmological constant, universe expands to infinite 
volume. With the massless scalar field playing the role of internal time, the infinite volume is 
reached in a finite value of internal time, <j> = <j)f. Its one consequence is that the quantum 
difference operator Q\ is not essentially self-adjoint. However, by choosing a self-adjoint extension, 
quantum evolution can be continued beyond 4>f. Extensive numerical simulations show that the 
detailed physics depends weakly on the choice of self-adjoint extension. The quantum difference 
operator 0a has discrete eigenvalues, however they depend on the choice of self-adjoint extension. 
Numerical evolution of states reveal that big bang singularity is resolved in the past evolution, and 
the evolution continues beyond (fif where a cycle of contraction starts. Bounces of the universe 
in different cycles occur when total energy density p 0.41ppi anc k. The quantum evolution 
from <j) = —oo to <j) = oo is thus composed of infinite cycles of expansion and contraction. 
Further, numerical simulations find an agreement between loop quantum evolution and effective 
dynamics. It is interesting to note that a non-trivial physical Hilbert space exists only when value 
of cosmological constant is smaller than a critical value of A = 3/7 2 A 2 [2]. Note that in the 
presence of massless scalar field, if A becomes equal to critical value, energy density would become 
greater than p max ~ Q-^-Pvunck- This observation is also tied to the von Neumann stability of the 
quantum difference equation, which we discuss in Sec. IV. 

Inflationary potential: Quantization of 4> 2 inflationary potential in spatially flat isotropic LQC 
is based on common features with that of the case of positive cosmological constant discussed 
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above |T0(] . In this model, scalar field 4> is n °t monotonic, however it can be treated as an internal 
clock locally. At each instant 4>, the quantum difference operator G is equivalent to the one in 
the case of positive cosmological constant, and is hence not essentially self-adjoint. For evolution, 
one thus has to carefully choose self-adjoint extensions between different 'time' slices. Preliminary 
numerical analysis with states which lead to a macroscopic universe at late times with this model 
indicates existence of bounce at p ~ 0.41/5pi anc k, and validity of the effective Hamiltonian approach 
till the scale of bounce. 

Bianchi-I spacetimes: Recently, a loop quantization of Bianchi-I model with a massless scalar 
field has been rigorously performed by Ashtekar and Wilson- Ewing (23|. Due to three spatial 
triads, the structure of non-singular quantum difference operator is much richer in comparison to 
the isotropic model. In this model, it is a partial difference equation in geometrical variables. Due to 
the underlying complexity of the partial difference equation, numerical evolution of physical states 
in this model are yet to be performed. One of the technical hurdles arises due to non-uniformity 
of discrete variables, which we discuss in Sec. IV. 

Numerical simulations have been performed for vacuum Bianchi-I model with a different quan- 



tization proposed earlie r by Chiou [24|]. Though this quantization suffers from fiducial cell scaling 
and infrared problems [57|, it nevertheless serves as a good playground to develop and test nu- 
merical methods beyond isotropic description. Detailed properties of the quantum theory, such 
as essential self-adjointness of the quantum difference operator, were found and techniques to ob- 



tain physical solution were introduced in Ref. 71J. In this analysis, one of the geometric variables 



plays the role of internal time, and relational dynamics between different geometric variables shows 
absence of singularity and existence of bounce. States remained sharply peaked throughout the 
evolution. It is to be noted that dynamics obtained from effective Hamiltonian constraint is shown 
to be a good approximation to the quantum evolution. 



IV. VON-NEUMANN STABILITY ANALYSIS OF DIFFERENCE EQUATIONS 

A useful way to understand various properties of the quantum difference operator 6 in LQC 
at large volumes, and its relation to the corresponding operator O in the Wheeler-DeWitt theory 
is by using ideas of von Neumann stability analysis. In various branches of physics, this analysis 
plays an important role in obtaining faithful numerical solutions of partial differential equations 
(PDE's) by discretization. Using von Neumann analysis, it is possible to determine whether a 
particular discretization of a PDE is stable, which by Lax-Richtmyer equivalence theorem is a 
necessary and sufficient condition for a consistent finite difference scheme to be convergent. Here 
consistency and convergence are used in the conventional sense of numerical methods (see for eg. 
[72I]). By consistency we imply that a smooth solution of a PDE is an approximate solution of 
finite difference equation obtained by its discretization when discreteness parameters go to zero, 
and convergence implies that a solution of finite difference equation approximates a solution of the 
PDE. 

In a conventional von Neumann analysis, for a given discretization of a PDE, say 1+1 dimension 
in space and time, one proceeds by a Fourier decomposition of the solution of finite difference 
equation in space which results in a recursion relation in time. Using this algebraic equation one 
can determine the amplification factor g between the values of the solution at two neighboring 
time steps. Von Neumann analysis requires that \g\ < 1 for all the roots of the algebraic equation. 
There are some important differences in the conventional usage of von Neumann stability analysis 
in various branches of physics and in LQC. These are noted below: 

1 . Unlike in the conventional applications of von Neumann analysis where differential equation 
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is fundamental, in LQC it is the difference equation which is fundamental. The discretization 
in geometric variable(s) is fixed by the underlying quantum geometry. Thus, in order to nu- 
merically solve the quantum Hamiltonian constraint in LQC, there is no freedom in changing 
the Planck scale discretization of quantum geometry to enhance numerical efficiency. 11 As 
emphasized in Sec. I, in the von-Neumann stability analysis in LQC, our goal is to compare 
solutions of difference equations which are uniformly discrete in different geometrical vari- 
ables obtained from loop quantizations with the Wheeler-DeWitt theory at large volumes to 
understand the infrared behavior of quantum difference equation in LQC. It is to be noted 
that as in conventional numerical methods, there is freedom to choose the discretization in 
(internal) time labelled by 4> appropriately. 12 

2. In the von Neumann analysis, stability issues are analyzed in the limit of discretization 
parameters going to zero. Since the discretization parameter(s) of geometry in LQC are 
fixed, one uses von Neumann analysis in the regime where the discreteness scale can be 
ignored. As an example, for the isotropic model with massless scalar, one can consider a 
large volume limit v ^> lp\ or v — > oo. In early works in LQC, different types of limits were 
used to analyze stability. Above criteria was first used in the work of Cartin and Khanna 



12] 



There is a difference on the roles of space and time in conventional usage of von Neumann 
analysis and in LQC. As discussed above, conventionally one performs a Fourier transform 
in space and analyzes the temporal behavior of amplification factor g. In isotropic LQC with 
massless scalar field as internal time, the Fourier transform is performed in </>, and spatial 
behavior of amplification factor at large volumes is analyzed. Even for other models so far 
studied in LQC, such as Bianchi spacetimes and Schwarzschild interior, one identifies one 
of the geometric variables as 'time', and analyzes the behavior of amplification factor in 



another geometric variable (see for eg. [74H7 1 



4. Strictly speaking, von Neumann analysis is applicable only for the case of constant coeffi- 
cients. For variable coefficients, as we will encounter in LQC, one can apply this method 
using frozen coefficient approximation. In conventional von Neumann analysis, this implies 
that one considers stability in a small neighborhood where coefficients in the PDE can be 
treated as constants up to a leading order. Extra care is needed to interpret results on 
stability in such a case. 

In LQC, the first work incorporating ideas of von Neumann analysis was by Bojowald and Date 



to understand properties of difference equations resulting from loop quantization 32J. At the stage 
of that early work, knowledge of the inner product and physical Hilbert space, which plays a vital 
role in ruling out unphysical solutions was not available. Rather, consistency of difference equation 
with respect to the Wheeler-DeWitt equation was analyzed, to rule out inconsistent quantizations. 
This analysis gave primary importance to the differential equation, as in the conventional usage 



1 Note that there are two types of discreteness studied here in LQC. Apart from the discreteness in quantum 
geometry resulting from the properties of area and volume operators in LQG, the other discreteness deals with 
structure of difference equation - whether it uniformly discrete in volume, area or any other geometric variable. 
(For a discussion on similar issues in LQG and broader context, we refer the reader to Ref. [z3]). A quantum 
difference equation, irrespective of whether it is uniformly discrete in volume or not, shares the same kinematical 
properties of discrete quantum geometry identified by the eigenvalues of the geometrical operators. 

12 This freedom is expected to be absent or severely limited if the scalar field is also polymer quantized using methods 
of LQG, instead of being Fock quantized. 
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of von Neumann analysis, and consistency of difference equation was verified in the limit of dis- 
cretization parameter becoming zero. This provided one notion of 'pre-classicality' in LQC. In a 
series of works, Cartin and Khanna investigated issues of 'pre-classicality' without taking the dis- 
creteness parameter to zero, but instead taking an appropriate limit to small spacetime curvature 



(or large values of volume) [12|, |35|, [38J. Solutions having pre-classical behavior have a desirable 
property expected from solutions of a consistent quantum theory with correct infrared behavior. 
Such solutions do not display unphysical small scale oscillations in regimes far from the Planck 
regime. 

Below, we focus on the von Neumann stability issues in the context of spatially flat isotropic 
models for improved quantization in LQC 0], and compare it with the old quantization in LQC 
This provides a flavor of the way this technique has been used in LQC to analyze 



consistency of quantum difference equations I75I477J ] . Our discussion has some parallels with the 



discussion of isotropic closed model in Ref. [59(, and on properties of eigenfunctions in presence 
of cosmological constant in Ref. [77J, where problems with infrared behavior have been discussed. 
After discussion of von Neumann stability issues for the spatially flat isotropic model, we summarize 
main results for other spacetimes in LQC. 

A. Stability analysis for flat isotropic model 

Results on the stability analysis for the improved quantization of spatially flat isotropic model 
with a massless scalar field considered in Sec. II were first reported by Cartin and Khanna [78]. 
Recently a comparison of improved and old quantizations in the isotropic model in perspective 
of behavior of solutions at large volumes has been made in the presence of cosmological constant 



771 ] . which can be understood using von Neumann stability. In the following we show the way 
stability analysis is performed in different models of LQC and contrast between improved and old 
quantization schemes. As we discuss below, important lessons and pitfalls of using this analysis can 
be understood by considering different matter. We first discuss the stability of quantum difference 
constraint as discussed in Sec. II for a massless scalar field and a positive cosmological constant. We 
then show that though the quantum difference equation for old quantization is stable for massless 
scalar, instabilities develop on inclusion of positive cosmological constant. As emphasized before, 
we note that in LQC it is the quantum difference equation which is fundamental, and stability or 
instability used below should not be confused with the conventional usage of terms in von Neumann 
analysis where PDE is fundamental and the goal is to determine whether a discretization of the 
PDE is stable or unstable. 



1. Isotropic model with massless scalar 

The loop quantization of massless scalar field in spatially flat isotropic model was discussed in 
Sec. II. Let us begin with rewriting of the quantum difference equation (I2.17P as follows, 

dp>{v,4>) = C+(i/)tf(i/ + 4A,<£) + C H<f (/;,</>) + C_(i/)#(i/-4A) (4.1) 

where coefficients C + , Co and C_ are given by 

C+{u) = ^v(y + 2\), C.( V ) = ^v(v-2\), C o (u) = -^0-u. (4.2) 

Our goal is to understand the way solution of the above equation changes (or amplifies) at neigh- 
boring steps at large volumes. To find this amplification factor, an equivalent short-cut to Fourier 
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transform procedure used in von Neumann analysis is to consider an ansatz (in the large volume 
limit) 

where n is an integer label of the discrete steps, and g is the amplification factor: (v + AX) / ^ (v) = 
g. Substitution of (|4.3p in (|4.7p yields the following quadratic equation in the amplification factor: 



C + {v)g 2 + (C (u) - u 2 )g + C-(y) = . (4.4) 

Since we are interested in the large volume limit, we consider the behavior of coefficients in a small 
neighborhood at large volume, where they are approximated as 

C + M^» = ^^O0 , (4.5) 

and 

Co-^c -|^, 2 + 0(1). (4.6) 

Von Neumann analysis requires both the roots of the quadratic equation (|4.4p to have modulus 
bounded by unity. Substituting (14.51) and (14. 6j) in eq. (|4.4l) . it is straightforward to see that the 
two roots for the amplitude g are equal to unity. Thus, by von Neumann analysis, the difference 
equation (|2.17p is stable at large volumes. The boundedness of the roots of the quadratic equation 
(|4.4p . suggests that there will be no sharp change in the solution at large scales. Recall that at large 
volumes, quantum difference equation for this model leads to the Wheeler-DeWitt equation (I2.2ip . 
Hence, according to von Neumann stability analysis, eigenfunctions of LQC at large volume should 
approximate those of Wheeler-DeWitt theory. This behavior is confirmed by Fig. 1 in Sec. Ill, 
where we find that eigenfunctions of LQC smoothly approach to those of Wheeler-DeWitt theory. 



2. Isotropic model: massless scalar with positive cosmological constant 

A rigorous quantization of positive cosmological constant in LQC has been performed recently 
by Pawlowski and Ashtekar, and extensive numerical studies of the resulting difference equation 
have been performed [3] (see also Ref. [3] for detailed properties of quantum difference operator). 
This case is interesting due to various physical and technical reasons, one of it being that since 
massless scalar has an equation of state —1, and that of cosmological constant is +1, the matter 
Hamiltonian spans the whole range of equation of state for matter satisfying weak energy condition. 
Thus, this model provides a formidable robustness test of the success of the quantization procedure 
in LQC. We now discuss the stability issue for this model on the lines of massless scalar field case 
discussed above. (A similar analysis of properties of solutions of quantum constraint operator with 
a different choice of lapse, has been performed by Tanaka et al [771] ). 

In the case of a positive cosmological constant with a massless scalar field, the quantum Hamil- 
tonian constraint takes the form of the following difference equation, 

8p{v, <P) = C+{y)V{y + 4A, <f>) + (C (u) + tt^GKu 2 )^^, 0) + C_ (i/)#(i/ - 4A) , (4.7) 

where C+ , Co and C_ are given by eq. (14. 2D , and are slightly different from the analysis in Ref. 
[l~9| due to factor ordering. Using the ansatz (|4.3j) in large volume limit, and following the steps 
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outlined in Sec. IVA1, we obtain a similar quadratic equation as (|4.4p . with the only change in the 
coefficient of the linear term in amplification factor: 

C Q (y) + WGKu 2 -u 2 ^ (l - 2 -^-) + 0(1), where A crit = -JL . (4.8) 

2X 2 V Acrit/ 7 * 



Using above equation and (|4.5]1 in the quadratic equation for g, we obtain the following roots: 

2 \ V2 

9= (1-2-r^) ±1 (1-2^1 -1 • (4-9) 




For A = 0, the roots of the equation are given by the massless scalar case discussed above. For 
A = Adit, the amplification factor gives two roots, equal to -1. For < A < Acrit j both roots are 
oscillatory with magnitude equal to unity. However, for A > A cr i t , one of the roots has magnitude 
greater than unity. Hence von Neumann analysis implies that for A > A cr it, eq. (j4.7j) can not be 
regarded stable at large volumes. 

This observation is interesting. As discussed in Sec. IIIC, it has been shown that for a loop 
quantization of massless scalar field with cosmological constant, a non-trivial physical Hilbert space 
exists only for A < A cr jt 0, Qjl]. The physical Hilbert space is zero dimensional for A > A cr j t . 
The slight discrepancy on the viability of the case A = A crr t between the analysis of Ref. 0,13 
and von Neumann analysis above, results due to the approximation involved in eq. (|4.8|) . which 
essentially ignores the role of massless scalar in the vanishing contribution of uj 2 term. If this 
order of contribution is included, von Neumann analysis leads to the same constraints on A as 



obtained by mathematical analysis of physical Hilbert space in Ref. [181. Il9l|. Further, the value of 
Acrit obtained from above stability analysis is the maximum allowed value of cosmological constant 
consistent with the maximum allowed energy density (p max = 0.41/?pi anc k) in massless scalar field 
case for a dense subspace of the physical Hilbert space in LQC. In the numerical simulations of 
semi-classical states for massless scalar with positive cosm olog ical constant, bounce of the scale 
factor has also been found to always occur at above value [l9j. It will be interesting to explore 
deeper relation between above observations originating from different directions in this model. 



3. Old quantization: massless scalar and cosmological constant 

We now discuss the von Neumann stability issue for the quantum difference equation which 
originates in the old quantization of LQC. In this quantization of LQC, the parameter p in eq. ([2.14|) 
was assumed to be a constant: p = p Q . With this assumption, the resulting quantum difference 
equation turns out to be uniform in eigenvalues of triad operator (which are proportional to area), 
rather than the volume operator jH]. A careful analysis of the eigenfunctions of thi s q uantization 



revealed problems with the infrared behavior for the case of cosmological constant [ill ]. Problems 
with the large volume behavior were also reflected in an early quantization of k = 1 model in LQC 
by Green and Unruh [lOj. As discussed in Sec. IIIC, a rigorous quantization with a massless scalar 
field in spatially flat model was performed in Ref. j3], which demonstrated a bounce. However, 
critical density turned out to be proportional to inverse of p^, and thus, was not invariant under 
rescaling of fiducial cell. Insights on the problems of this quantization using effective dynamics are 
discussed in Sec. V. 

In the following we first do analysis with a massless scalar, and then include a positive cos- 
mological constant. Unlike the case of the improved quantization in isotropic LQC, results in the 
presence and absence of positive cosmological constant turn out to be starkly different. (For this 
analysis, we consider lapse function to be unity, and we work in the factor ordering of Ref. |7y with 
the same value of discreteness parameter \i a = Sy/3/2). 



21 



In the presence of massless scalar field, loop quantization with \x = \i yield the following 
quantum difference equation [3] 

8ttGB(p) djHM) = (/+(/*) *(0 + 4a*o,$ + fo(n)*M) + /-(M)*(M-4Mo)) (4.10) 
where 

/+(/*) = 2 V "6~ (7//) 3 / 2 ' /i + 3/x °| 3/2 ~ I A* H- Mo 1 3/2 , /-O) = /+(//-4/i ) and / (/x) = -/_(//) 

™° (4.H) 

and -B(p) denotes eigenvalue of inverse volume operator l/|p| 3 / 2 , given by 

B ^(4) I ' ! (s) f (" + ^ H ^ ft|1 ' , ) < (4 - l2) 

We consider the regime /i very large compared to fi a , and using the ansatz ^(n,<p) = g n l^o^4> j n 
von Neumann analysis, we obtain the quadratic equation 

Ui^g 2 + (/o(//) - ^GB{p)h 2 u 2 )g + /_» = . (4.13) 

The coefficients of this quadratic equation can be approximated for large \x as, 

3 l Pl 1/2 1/2 



- /-(") = 2 VT^I^ + } (414) 



and 



/o - 8vrGB(p)u, 2 ~ ^M^h ^ + 0(^ 1/2 ) • (4.15) 

Substituting (|4.14p and (|4.15|) in eq. (|4.13p . one finds that both the roots of the difference equation 
(|4.16|) are equal to unity. Thus, for the case of massless scalar field, old quantization appears to 
yield a similar result for stability analysis at large volumes as the improved quantization. Since the 
domain of validity of von Neumann analysis is only for the large volumes, it will be a mistake to 
extend this conclusion to small volumes. As we noted above, and will further discuss in the next 
section, though the infrared behavior as analyzed from the perspective of von Neumann analysis is 
similar for two quantizations, old quantization in LQC for massless scalar field leads to unphysical 
quantum gravitational effects which are easily ruled out. 

We now turn to the interesting case of cosmological constant in the old quantization of LQC. 
Since this case is known to be problematic in the large volume regime [ll|, we expect von Neumann 
analysis to yield amplification factor with modulus greater than unity. In the presence of positive 
cosmological constant, the quantum difference equation turns out to be 



^GB{p)dl^{^4>) = /+(m)*O*+4Mo,0)+ /o(m) + 2 



(4.16) 

where /+, /o and /_ are given by eqs. (|4.1ip . As before, the von Neumann stability analysis yields 
a quadratic equation for amplification factor, which is 



U^)9 2 + ( fo ~ 8nGB(p) hW + 2 (^^j^A^) g + f_M = 0. (4.17) 
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In the regime where \x is large, /+ and /_ are approximated by eq. (|4.14j ). The coefficient of the 
linear term of the above quadratic equation approximates to, 

fo ~ 8nGB(p)^ + 2 ^ ^ - 2A ^ ^ + 0{pT^) . (4.18) 

Using these expressions in (|4.17|) . within the order of approximation of large fi, we obtain following 
two roots 

51 = 0, & = -^Allfrno) 3 ti . (4.19) 

Since 7 and fi are fixed, for any given value of A, there exists a sufficiently large value fj,, such that 
I c/2 1 > 1- Thus, this model fails the von Neumann stability test at large volumes. This exercise 
gives a valuable lesson that one must be careful in claiming stability of a quantization scheme, and 
it is important to understand behavior of the difference equation with different matter. The above 
limitation of old quantization has been understood in literature in various ways, and is reflection 
of the bad infrared behavior of the quantization, and is thus ruled out. 

We conclude this subsection with two remarks, one_on problems with arbitrary discretizations 
which are sometimes proposed in loop quantization and other on the closed models. 

Remark 1: In literature, arbitrary discretizations have also been considered in lattice refined 
models inspired by LQC. 13 For the isotropic models, these are based on considering holonomies of a 
general variable introduced in Sec. IIA, P m = c\p\ m with —1/2 < m < It is straightforward 
to repeat the above analysis for any value of m in the above range for cosmological constant, 
and one encounters similar difficulties as noted above for the old quantization in LQC. In fact, 
in isotropic LQC, it is only for the case of m = —1/2, which corresponds to P m = b, that the 
von- Neumann stability analysis gives positive result for the case of positive cosmological constant. 
Thus, the discretization in quantum theory is highly restrictive. This observation is in agreement 
with the one on uniqueness of discretization in LQC obtained in Ref. [13] using effective dynamics 
(discussed in Sec. V). Recall that in Sec. IIA, we discussed that for m 7^ —1/2, one does not 
obtain invariance of phase space variable P m . We will discuss in the next section, how these 
different observations tie up with each other. 

Remark 2: As mentioned earlier, based on an early quantization of isotropic models in LQC, 
Green and Unruh investigated properties of quantum difference equation and the behavior of 
eigenfunctions in spatially closed k = 1 model with a massless scalar field fioj] . In this model, 
classical GR predicts a recollapse of the universe at a volume determined by the initial energy 
density of the scalar field. If the quantum theory leads to a consistent classical behavior at small 
spacetime curvature, one expects wavefunctions to exponentially decay near the volume at which 
recollapse occurs. Authors of the above work, did not find this behavior. Instead, they found that 
in general, solutions grow exponentially at large volumes. It is to be noted that at that stage of 
analysis, a rigorous quantum theory of spatially closed model in LQC was lacking. In particular, 
little was known about the inner product which could provide insights on physical and unphysical 
solutions. Nevertheless, Green and Unruh's work brought to light limitations with the earlier 
quantization in LQC. As discussed in Sec. Ill, these challenges were overcome successfully in an 



On similar lines, ad-hoc modifications to the Hamiltonian constraint and their effect on bounce have also been 
investigated 79]. Though such studies may be useful to understand the way properties of difference equations 
are affected by certain terms, being not derived from loop quantization, such studies give little insights on the 
robustness of bounce in LQC. 
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improved quantization which demonstrated correct ultraviolet and infrared behavior 14| . Extensive 
numerical simulations showed that a spatially closed LQC universe, does indeed recollapse at the 
volume predicted by the classical theory, and resolves the singularity via a non-singular bounce. 
From the perspective of discretization, it is to be noted that the quantization Green and Unruh 
used is uniformly discrete in triad eigenvalues, whereas the correct quantization leads to a difference 
equation which is uniformly discrete in volume. Recently, a von Neumann stability analysis of 
quantization considered by Green and Unruh shows that one of amplification factors has modulus 
greater than unity [o^ ]. In contrast, both of the roots in quadratic equation for amplification 
factors in the quantization of Ref. [1J], have modulus equal to unity, as in the flat isotropic model. 



B. Anisotropic spacetimes and Schwarzschild interior 



In this part we the summarize results on understanding of properties of quantum difference equa- 
tion, including von Neumann stability issues for loop quantization of spacetimes beyond isotropy 
assumption. As discussed earlier, in the early stages of investigations in LQC, a rigorous formu- 
lation of quantum theory, including the physical inner product, was not available, and resulting 
physics was not well understood. In absence of these tools, various early works used notions of pre- 
classicality at large volumes to understand properties of quantum difference equation in relation to 
GR. Pre-classicality was obtained in different ways, depending on the procedure of taking underly- 
ing limit, whether of discreteness parameter [32], or to the small spacetime curvature 1_2|, ESS. 
In the above approach, usage of gen erating function technique was particularly useful to obtain 
pre-classical solutions analytically [33t] . With generating functions method, Cartin and Khanna [l~2| 
showed non-existence of pre-classical solutions for an early quantization of vacuum Bianchi-I model 
proposed by Bojowald [4|. It was then argued by Date that if one does not require pre-classical 
solution to correspond to the sector with a vanishing volume, then pre-classical solutions can exist 
[3^ |. Cartin and Khanna then studied this model using separation of variables and investigated 
existence of pre-classical solution for all values of triads. It was found that existence of such pre- 
classical solutions depends on the sign of the separation parameter [IH]. It should be noted that 
presence of matter or cosmological constant may alleviate the problem of pre-classical solutions, as 
was shown in the analysis of Ref. 80(] . For the same quantization in the case of vacuum Bianchi-I 
LRS model a von Neumann stability analysis with a symmetric as well a non-symmetric forms of 
the quantum constraint revealed instabilities in large regions of spacetime 74]]. In contrast, the 
quantum difference equation for an earlier quantization of Bianchi-I vacuum spacetime [24j] , written 
in the separable form turns out to be stable [8l|. Generating functions were also used to study 
pre-classical solutions in the Schwarzschild interior and limitations of the quantization proposed 
in Ref. |27j were identified 38] . A von Neumann stability analysis of the same quantization of 



Schwarzschild interior shows instabilities [74]], whereas a different quantization of Schwarzschild in- 
terior heuristically motivated by Boehmer and Vandersloot [2^] demonstrates stability 75]. Nelson 
and Sakellariadou, studied von-Neumann stability of the quantum difference equation for Bianchi-I 
vacuum spacetime resulting from the analysis of Ashtekar and Wilson- Ewing [23|] , and claimed that 
in the explicit implementation, the difference equation is unstable [76(. As we have discussed in 
Sec. IVA and above, it can be naive to extend conclusions from one particular choice of matter 
(or absence of it) to general matter. It is quite possible that in presence of massless scalar the 
quantum difference equation turns out to be stable [82|. Since some of the detailed properties of 
the quantum theory, such as self-adjointness of quantum constraint operator and the associated 
inner product, in above Bianchi-I model remain to be fully understood, it is also possible that 
instabilities correspond to unphysical solutions which are weeded out in the physical Hilbert space. 
The effect of implementation of quantum difference equation in an implicit scheme, or by use of 
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different variables remains to be understood. 

We now discuss an important technical issue which arises in the quantization of Bianchi space- 
times and the Kantowski-Sachs model. In these models quantum constraint operator in LQC is 
a difference equation in more than one variables. Instead of an ordinary difference equation as in 
the case of isotropic models, one has to solve a partial difference equation. This partial difference 
equation can not be expressed in a form such that it is uniformly discrete in all the variables. Using 
appropriate redefinition of variables, it has been shown that the partial difference equation can be 
made uniformly discrete in one of the variables (see for eg. Ref. [75], and Ref. [2^] for a rigorous 
implementation at the quantum theory level in Bianchi-I model). Due to the non- uniform step 
size, the structure of difference equation is such that the value of the wavefunction at some future 
discrete step can not be recursively obtained from previous iterations. Further, the discretization 
step in one direction can also depend on the one in the other direction, and the partial difference 
equation is not variable separable. 

To overcome these problems, Sabharwal and Khanna have recently used a local interpolation 
method J39(], and Nelson and Sakellariadou have proposed a Taylor expansion scheme for interpo- 
lation [40(]. We illustrate below the difficulty of non- uniform grids and summarize the interpolation 
method of Ref. [39] in a simple case where separability in variables is allowed. Let us assume a 
quantum constraint with two discrete variables (<5, r), which are separable, and the partial differ- 
ence equation can be written as two coupled difference equations. If Cj (where i = 0, 1, 2) denote 
coefficients in the r difference equation, then say at 5 = S , recursive relation can be used to 
determine C2, given C Q and C\ as initial data. However, due to non-uniform stepping, at 5 = Si, 
the recursive relation in r would have coefficients which depend on value of r not obtained from 
uniform stepping of the values at 5 = 5 , such as C2-E, Ci and C2+ t - Using difference equations, it 
is not possible to determine C2- t from C2, C\ and C , and hence one can not use recursive relation 
at 5 = Si to obtain C^+e- In Sabharwal and Khanna's algorithm, one performs a least square fit at 
S = S to obtain a locally accurate formula which is used to determine C^-e from Ci,C\ and C Q . 
This allows to compute C2+ e at 5 = Si using difference equation at S = Si. This value is used to 
obtain a revised locally accurate formula in the neighborhood of r = 2 — e, 2, 2 + e, which is then 
used to compute coefficients in the next step S = 5%- The local interpolation method is applicable 
in the regime where solutions are varying slowly in the neighborhood. Due to similar behavior of 
solutions at small spacetime curvature (or large volumes) in various models, above method is useful 
to understand large volume behavior of quantum difference equations, von Neumann stability and 
evolution of semi-classical states at late times. However, in the regime where solutions can change 
rapidly in the neighboring points, a more general method, such as the one based on Taylor series 
expansion for interpolation on neighboring points is expected to be more reliable, as it requires 
only analyticity of the solution and provides a better control on the accuracy of method [40]. At 
this stage, more work is needed to understand the domain of applicability of these methods. 

V. EFFECTIVE DYNAMICS 

In the previous sections we have discussed different numerical methods to understand detailed 
features of quantum difference equation in LQC. Interestingly, various aspects of quantum dynam- 
ics in LQC can be understood also in terms of an effective spacetime description using effective 
dynamical differential equations derived from an effective Hamiltonian constraint. This description 
is derived in LQC using geometric formulation of quantum mechanics where one treats the space 
of quantum states as an infinite dimensional quantum phase space with a symplectic structure and 
relates it to the classical phase space by seeking a faithful embedding [83]. Such an embedding can 
be obtained by a judicious choice of states, such as coherent states which at late times peak on 
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FIG. 4: This plot demonstrates the validity of effective dynamical trajectory for the massless scalar model in 
LQC for the numerical simulation of quantum evolution of initial state discussed in Fig. 2. The expectation 
values of volume observable are peaked on the effective trajectory (solid line) obtain from eq . (|5.1|) throughout 
the evolution. 



classical trajectories in a macroscopic universe. In LQC, an effective Hamiltonian up to quantum 
corrections corresponding to properties of the states such as fluctuations, has so far been achieved 
for the spatially flat model in presence of a massless scalar [42J and dust [41] ]. and generalization 



for arbitrary matter is in progress 841 ] . Interestingly, using extensive numerical simulations on 



evolution of states with quantum difference equations in LQC, validity of effective Hamiltonian 
has been verified for models in the presence of spatial curvature 14|, ll_6|], inflationary potential 



701 ] and even anisotropics [Til . |85[|. These works, hence, strongly indicate that the validity of 
effective Hamiltonian approach in LQC is more general. We illustrate the comparison of effective 
dynamics with quantum evolution in LQC for the numerical simulation discussed in Fig. 2, in Fig. 
4. It can be clearly seen that effective dynamical trajectory obtained from the effective Hamilto- 
nian constraint discussed below (eq. (l5.ip ) is in an excellent agreement with the expectation values 
of volume observable in LQC computed by state evolution, at all scales. In a recent work, this 
agreement has been observed even for states which bounce at very small volume [66]. 

In LQC, effective dynamics has yielded many insights on the new physics at Planck scale, in- 
cluding issues ranging from genericity of singularity resolution J54l456| , and phenomenological issues 
such as attractor behavior in inflationary dynamics 
of inflation 



43l . |45| . 1471 ] and implications for probability 



45l . |46|], to detailed analysis of approach to Planck scales in anisotropic models in 



LQC 5ll. |52[|. and on inclusion of inhomogeneous degrees of freedom [48] (for discussion of various 
applications see, Sec. V of Ref. [1]). Various results obtained in effective dynamics rely on a 
use of numerical methods to solve coupled ODE's, which carefully incorporate control on error 



Another approach which can be used to obtain an effective Hamiltonian is a truncation method [86j]. This method 
is also based on geometric formulation of quantum mechanics, and is similar to the approach of order by order 
perturbative analysis. It allows a discussion of broader set of state, however extra care is needed to systematically 
control truncation errors. This method has not been as widely tested in LQC models as the one based on embedding 
approach. For more discussion of analysis based on this approach we refer the reader to Refs. [87l . [j^ . 
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terms, arising for example in a slight non-vanishing of the Hamiltonian constraint in a numerical 
simulation. For examples of works with detailed discussion on these issues in LQC, we refer the 



reader to Refs. [45|, |47J, [5l|, |K 

Amongst various applications of effective dynamics, a very useful one concerns with its use to 
gain insights on the problems with certain discretizations in LQC, such as old quantization 0, 0| 
or arbitrary discretization as proposed in lattice refinement scheme [62]. Below we discuss this 
aspect of effective dynamics, and show the way it leads to reveal similar problems with certain 
quantizations as demonstrated in von Neumann stability analysis. Our discussion is based on 
the analysis in Ref. [13]. Since this part of the manuscript has some parallel with the classical 
cosmology, we work with the lapse N = 1. The effective Hamiltonian constraint for a spatially flat 
isotropic model with a massless scalar field in LQC is given by [42] 

3/i sin 2 (A6) „ , „ ,~ „ s 

--^-Lv + HtnO, (5.1) 

where = p^/2\p\ s ' 2 is the Hamiltonian for the massless scalar. It is straightforward to verify 
using eq. (|2.4p . that in the regime Xb <C 1 above equation yields the classical Hamiltonian constraint 
(|2.8p for lapse N = 1. Recall that on classical solutions of GR for massless scalar model, b = oc 
o -3 , and hence, effective Hamiltonian constraint approximates classical GR at large volumes. The 
same conclusion is reached for effective Hamiltonian constraint with arbitrary matter satisfying 
weak energy condition [l^. This shows that effective Hamiltonian constraint (|5.1|) leads to a 
consistent infrared behavior. To understand the ultraviolet behavior, we consider the expression 
for energy density of the scalar field, which using eq. H5.ip . turns out to be, 



IP* 3 sin 2 (Afr) 
P ~2V^~8ttG 7 2 A 2 ■ 15 ] 

The maximum value of energy density is given by p max = 3/87rG7 2 A 2 0.41^^^, a universal 
constant. Recall that this is the same value where bounce occurs in numerical evolution of quantum 
states in LQC 0] and the value which emerges as an upper bound on expectation values for all 
states in the physical Hilbert [l| . It turns out that the effective dynamics also predicts a bounce at 
this energy density, and a period of superinflation between p max < p < p ma , x /2 [gij, whose many 
interesting implications have been studied in LQC. 

So far we have discussed the effective Hamiltonian constraint in the improved quantization in 
LQC which yields a quantum difference equation which is uniformly discrete in volume. We now 
discuss some features of the effective dynamics in cases when quantization leads to a difference 
equation not uniformly discrete in volume. Examples are the cases of old quantization in LQC 
where the quantum difference equation is uniformly discrete in triad eigenvalues [tJ, and lattice 



refinement scheme where arbitrary discretizations are allowed [62|. All these cases can be studied 
by considering phase space variable P m — c\p\ m , for -1/2 < m < (introduced in Sec. IIA), and 
its conjugate x = |p| 1_m /(l — rn). For a quantization based on these variables, the expression for 
energy density from the effective Hamiltonian constraint, which is similar to (|5.ip with sinusoidal 
term given by sin 2 (Ap fl P 9 ), turns out to be fl~3| 



8vrG 



(2m+l)/(2m-2) 



where Xp m is a constant related to the area gap in LQC. Note that except for the case of m — —1/2 
which corresponds to P m = b (and a uniform discretization in volume), the maximum value of p 
depends on p^. As discussed in Sec. IIIC, this leads to a serious problem because ps is not invariant 
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under the rescaling of the fiducial cell, and the maximum energy density in such quantizations, at 
which bounce occurs, can be reached at arbitrarily small spacetime curvature by a simple rescaling 
of the fiducial cell. Unfortunately, even fixing a compact spatial topology such as a 3-torus, with 
say V = 1, in these quantizations does not help and problem for a consistent infrared behavior 
persists. Note that for arbitrary m, P m varies with scale factor on classical trajectories as, 

P m = cp m = 7 a a 2m oc a -(3»+i-4m)/2 ( 5 4 ) 

where we have used classical equation d oc a^~ 3w+1 ^ 2 following from classical Friedmann equation 
and stress-energy conservation law. For — 1/2 < m < 0, P m increases with an increase in the scale 
factor of the universe for — 1 < w < 0. This implies that the term sin 2 (Xp g P g ) in the effective 
Hamiltonian constraint for — 1/2 < m < leads to large deviations from GR even for matter which 
satisfies weak energy condition. This result is in confirmation with our discussion on problems with 
old quantization (which corresponds to m = 0)and arbitrary discretizations (—1/2 < m < 0) for 
von Neumann stability issues of quantum difference equation. In the light of the above argument, 
in Sec. IVA3 one would have expected that the old quantization of LQC will be problematic in 
achieving correct classical limit for the case of cosmological constant, which turns out to be true. 
In fact, it is only for m = —1/2 that one obtains a consistent infrared behavior. Hence, out of all 
possible discretizations, uniform discretization in volume is unique in this sense [liil]. 15 

Above discussion shows that effective Hamiltonian approach can be very useful in predicting 
detailed features of quantum difference operators which can be confirmed by more sophisticated 
numerical methods, such as von Neumann analysis. In particular, innocuous rescaling of variable 
under change in fiducial cell can teach us important things. In the same spirit as above, different 



quantizations in Bianchi-I models have been discussed 57(, which show the viability of the recent 
quantization by Ashtekar and Wilson-Ewing j^] in comparison to an earlier quantization by Chiou 
24]. It will be interesting to see if the analysis of Ref. [57( provides useful insights on the von 



Neumann stability analysis of Bianchi-I models in LQC with matter. 



VI. SUMMARY AND OUTLOOK 

In this manuscript we have reviewed various numerical techniques which have been used in 
LQC to understand the behavior of quantum difference operator and physical states. Since the 
applied numerical methods are diverse, we used the quantization of massless scalar field in isotropic 
spacetime as a platform to discuss them in a coherent way. After providing a brief summary of 
the classical framework and basics of quantization and how it leads to the quantum Hamiltonian 
constraint as a difference equation with a uniform discretization in volume, we described the 
underlying procedure of obtaining physical states using an FFT and with an internal time evolution. 
These techniques have been extensively used in various isotropic models and also in preliminary 
investigations of anisotropic models in LQC. In all the models which have been quantized so far in 
LQC, these methods reveal existence of bounce when the spacetime curvature becomes Planckian. 
A lot of work in LQC has been devoted to understand the behavior of quantum difference operator 
at large volumes and consistency of quantization at infrared scales by comparison with the Wheeler- 
DeWitt theoru using von Neumann stability analysis. Illustrating these techniques for the isotropic 
model, we showed the way improved quantization of Ref. [8j overcomes problems associated with 
the old quantization in LQC Insights on these issues were also discussed using effective 

dynamics, and we highlighted the complementarity of these approaches to understand infrared 



15 A different analysis based on factor ordering ambiguities results in a similar conclusion [Hit ]. 
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issues. It is notable that the unique viable discretization which is established from these works, 
emerges naturally from loop quantization §]. It will be useful to apply these ideas in greater 
synergy for anisotropic and black hole spacetimes to gain insights on the properties of quantum 
difference equations in these models. 

Numerical techniques in LQC have so far played a very important role in its progress, however, 
the field of numerical loop quantum cosmology is still in its infancy. Various avenues remain to 
be explored in near future which open opportunities for contributions from experts in numerical 
and computational methods. We have discussed some of the open directions in the manuscript, 
which include development of techniques to solve partial difference equations with non-uniform 
discretization carrying forward the works of Refs. 39, |4(j. These methods will be in particular 



helpful to uncover detailed physics of anisotropic and black hole spacetimes in LQC. A related arena 
is to explore the quantum properties of spacetime in the bounce regime, on the lines of numerical 
techniques developed in classical GR to explore the structure of spacetime on approach to space-like 
singularities 9CJ, |9l| . Work on these lines is expected to provide deep insights on the way quantum 
gravity affects mixmaster behavior, and critical phenomena in gravitational collapse. For the 
Bianchi models, these avenues are being investigated using effective dynamics which already provide 
an interesting result, that quantum geometric effects lead to Kasner transitions across the bounce 
in Bianchi-I spacetimes [52J. Critical phenomena and Choptuik scaling in gravitational collapse 
scenarios is also being studied but so far only partial effects from LQG have been incorporated (see 



for eg. [9214941] ). These works indicate that underlying quantum geometry may lead to a mass gap 



in gravitational collapse scenarios Whether or not these results are borne out in a rigorous 

quantization can only be answered by carrying out extensive numerical simulations in near future. 
Investigations on these lines are also expected to provide us with insights on the phenomena of 
black hole evaporation. Recently, insights from singularity resolution in LQC were fruitfully used 
to gain understanding of this issue for Callan-Giddings-Harvey-Strominger (CGHS) black holes 
[hH], and extensive numerical analysis revealed new insights on the quantum evaporation of these 
black holes [9^]. Future work on similar lines with Schwarzschild black holes is expected to be a 
fertile ground for research where numerical methods will play an important role. 

Going beyond the isotropic models, or asking more refined questions even in the isotropic models 
will require adapting algorithms to high performance computation techniques. To illustrate this, 
let us consider computational costs on an interesting avenue of research - to test the validity of 
effective dynamics in extreme regimes (i.e. at volumes comparable to Planck volume) with states 
which are more general than the ones considered in this manuscript. In the spatially flat isotropic 
model, note that the characteristic speeds in Wheeler-DeWitt equation (in the v representation), 
to which quantum difference equation in LQC approximates to at large volumes, are proportional 
to volume. Since the discreteness in volume is fixed by the quantum geometry, this implies that 
the 'time' step to numerically integrate equation in <j> determined by CFL condition is inversely 
proportional to the volume. Due to this reason, the computational cost of numerically evolving 
states with wider spreads is much higher with the same techniques as used for sharply peaked 
states. Note that a wider state, would require the boundary in volume to be considered at a much 
higher value than the state which is sharply peaked. And, an increase in boundary volume, implies 
choosing a smaller time step. As an example, for the numerical simulation used in generation of 
Fig. 2, the outer boundary was taken at v = 400000 and the time interval was divided in 300000 
steps. The simulation took approximately 45 minutes on a dual core Opteron 280 machine. A 
naive implementation of the same algorithm for an outer boundary, say at v = 4 x 10 10 (which is 
typically required for certain types of quantum states), increases the computational time by 10 10 
hours! These costs would only grow for asking similar questions in anisotropic models. Certainly, 
better algorithms harnessing the power of high performance computing are required, which are 
being developed [66] . These algorithms will provide us with tools to understand robustness of new 
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physics in the Planck regime and answer fundamental questions on the emergence of continuum 
spacetime from quantum geometry. Development of these algorithms will also open new avenues to 
perform extensive numerical evolution of states in anisotropic and black hole models using HPC's 
in the near future. 
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